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ABSTRACT 

The linear growth rate is commonly defined through a simple deterministic relation 
between the velocity divergence and the matter overdensity in the linear regime. We 
introduce a formalism that extends this to a nonlinear, stochastic relation between 9 = 
V • v(x,t)/aH and 6. This provides a new phenomenological approach that examines 
the conditional mean (0|5), together with the fluctuations of 6 around this mean. 
We measure these stochastic components using N-body simulations and find they are 
non-negative and increase with decreasing scale from ~10% at k < 0.2hMpc -1 to 
25% at k ~ 0.45hMpc -1 at z = 0. Both the stochastic relation and nonlinearity are 
more pronounced for halos, M ^ 5 x 10 12 Mgh _1 , compared to the dark matter at 
z = 0 and 1. Nonlinear growth effects manifest themselves as a rotation of the mean 
(0|5) away from the linear theory prediction — /lt< 5, where /lt is the linear growth 
rate. This rotation increases with wavenumber, k, and we show that it can be well- 
described by second order Lagrangian perturbation theory (2LPT) for k < O.lhMpc . 
The stochasticity in the 6-6 relation is not so simply described by 2LPT, and we 
discuss its impact on measurements of /lt from two point statistics in redshift space. 
Given that the relationship between 6 and 9 is stochastic and nonlinear, this will have 
implications for the interpretation and precision of /lt extracted using models which 
assume a linear, deterministic expression. 

Key words: Methods: N-body simulations - Cosmology: theory - large-scale structure 
of the Universe 


1 INTRODUCTION 

The clustering of galaxies on Mpc scales in the Universe is 
a fundamental cosmological observable which allows us to 
constrain key parameters of the ACDM model and to look 
for deviations from this standard model. Understanding the 
relationship between peculiar velocity flows and the large 
scale mass distribution is crucial to interpreting the cluster¬ 
ing signal measured in redshift space, where these velocities 
distort the clustering amplitude along the line of sight (see 
e.g. Peacock et al. 2001; Guzzo et al. 2008; Blake et al. 2011; 
Reid et al. 2012; Beutler et al. 2014). In this paper we in¬ 
vestigate the assumptions of a linear and deterministic rela¬ 
tion between the peculiar velocity and overdensity fields at 
a range of scales and redshifts. We present a general formal¬ 
ism where deviations from linearity and determinism can be 
viewed separately in the two point clustering statistics of 
the velocity divergence auto and cross power spectra. This 
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approach represents a new phenomenological tool based on 
a stochastic description of nonlinear effects. 

One of the key aims of future galaxy redshift surveys 
(Cimatti et al. 2009; Spergel et al. 2013; Eisenstein & DESI 
Collaboration 2015) is to measure this linear perturbation 
theory relation between the density and velocity fields, re¬ 
ferred to as the linear growth rate, to less than 1% preci¬ 
sion using the redshift space clustering statistics of different 
galaxy tracers. This level of accuracy has motivated a lot of 
work in developing a precise model for the two point cluster¬ 
ing statistics either as the correlation function in configura¬ 
tion space (e.g. Reid & White 2011; Bianchi et al. 2015) or 
the power spectrum in Fourier space (e.g Peacock & Dodds 
1994; Scoccimarro 2004; Jennings et al. 2011; Seljak & Mc¬ 
Donald 2011; Taruya et al. 2013). Note that many of these 
studies are based on a mix of assumptions of either a linear 
and/or deterministic density velocity relation. 

Current models for the two point clustering statistics in 
redshift space that include perturbation theory expansions 
have been shown to be an improvement over linear theory in 
modelling these redshift space clustering statistics. Although 
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all are limited to very large scales k < 0.15frMpc -1 at low 
redshifts (see e.g Scoccimarro 2004; Jennings et al. 2011; 
Kwan et al. 2012) and moreover may only apply to highly 
biased tracers (Reid & White 2011); none of the models can 
recover the linear growth rate to a percent level accuracy on 
the scales which will be probed by future galaxy surveys. If 
we are to limit our analysis of redshift space distortions to 
large scales, where quasi-linear theory models apply, then 
it is worthwhile investigating both where the assumptions 
of a linear and deterministic relation between the density 
and velocity fields breaks down and how well perturbation 
theory expansions can recover these components. 

This formalism involving the decomposition of the two 
point statistics into nonlinear and stochastic components is 
both well defined and consistent with a full perturbation 
theory expansion of all the nonlinear effects. The approach 
provides an alternative, more phenomenological description 
of such nonlinear effects. In considering either the galaxy - 
dark matter overdensity relation or the velocity- overden¬ 
sity relation, there is a general notion of stochasticity which 
is often not well defined and or vaguely explained as due 
to a nonlinear coupling of modes. In this paper, our use of 
the term stochasticity refers to the break-down of a deter¬ 
ministic relation that exists in the linear regime between 
the growing overdensity field and the velocity divergence. 
We also discuss the connection between such a notion of 
stochasticity and mode coupling in standard perturbation 
theory. 

It is well known that the halo or galaxy overdensity field 
does not trace the dark matter field faithfully and that the 
relation between the two is generally described by a linear 
bias term which is scale independent and is different for dif¬ 
ferent galaxy tracers (see e.g Dekel & Lahav 1999; Kravtsov 
& Klypin 1999). Recently there is renewed interest in consid¬ 
ering the stochasticity in this relation on large scales (Seljak 
& Warren 2004; Bonoli & Pen 2009; Sato & Matsubara 2013) 
where previously we would have assumed a linear, determin¬ 
istic relationship to hold. Also, as noted in Seljak & Warren 
(2004), dominant perturbative corrections come from mode 
coupling at wavelengths close to the wavelength of the mode 
itself. Long wavelength modes sampled from a finite volume 
can have significant fluctuations which would give rise to 
significant fluctuations in second order corrections. 

There have been many studies that have compared 
the two point statistics of the matter and velocity diver¬ 
gence fields and found them to be nonlinear on large scales 
(k ~ 0.1/i/Mpc) which are traditionally considered the lin¬ 
ear regime (Scoccimarro 2004; Percival & White 2009; Jen¬ 
nings et al. 2011; Jennings 2012; Crocce et al. 2012; Jen¬ 
nings et al. 2012). Jennings (2012) measured this nonlinear¬ 
ity as the deviation of the velocity divergence power spectra 
Pee '■= (0(k)0*(k')) and Pes : = (0(k)<5*(k')) from linear per¬ 
turbation theory predictions and found it to be at the level 
of 20% and 10% respectively at k ~ O.lhMpc -1 . Note that 
these nonlinear features are at the level of the ensemble av¬ 
eraged two point statistics. In contrast, in this work we wifi 
examine the velocity divergence - overdensity relation, 6 — 8, 
in Fourier space at each wavenumber where we can separate 
the notion of nonlinear and stochastic effects. 

Bernardeau et af. (1999) investigated the statistical re¬ 
lation between the density and velocity fields in the mildly 
nonlinear regime, focusing on the conditional probability 


distribution P(6\8) of the smoothed fields in configuration 
space. This study of the stochastic relationship between the 
two fields used simulations of a small volume, (200Mpc//i) 3 , 
and low resolution, 128 3 particles, by today’s standards. 
Given the high resolution and large volume simulations we 
have available today and our knowledge of how sensitive 
velocity statistics are to resolution effects (Pueblas & Scoc¬ 
cimarro 2009; Jennings et af. 2011, 2015; Zheng et al. 2014; 
Biagetti et al. 2014) it is important to revisit this study. In 
this paper we explore a formalism that defines both a nonlin¬ 
ear and a stochastic relation between the velocity divergence 
and the conditional mean value of this function at a given 
overdensity. We also investigate the variance of the velocity 
divergence around this relation as a function of scale, which 
defines a stochastic description of nonlinear effects. 

The paper is laid out as follows: In Section 2 we de¬ 
scribe the N-body simulations and tessellation techniques 
used to measure both the density and velocity divergence 
fields of dark matter and halos in this paper. In Section 3.1 
we present the linear perturbation theory relation between 
the density and velocity fields. In Sections 3.2 we outline the 
main formalism in this paper which defines the nonlinearity 
and the stochastic relation between the velocity divergence 
and overdensity fields and give expressions for the two point 
statistics. In Section 4 we present our results. The measure¬ 
ment of the conditional mean relation and scatter about this 
mean are presented in Sections 4.1 and 4.2 for dark matter 
and in Section 4.3 for halos. In Section 4.4 we relate the two 
point functions in this paper to both one loop standard per¬ 
turbation and second order Lagangian perturbation theory 
predictions. In Section 5 we discuss the impact of a nonlin¬ 
ear and stochastic relation between the velocity and density 
fields on models for the power spectrum in redshift space. 
In Section 6 we summarize our results. 


2 DENSITY AND VELOCITY TWO POINT 
STATISTICS FROM N-BODY SIMULATIONS 

In section 2.1 we present the details of the dark matter N- 
body simulations and the MultiDark halo catalogue used in 
this work. In Section 2.2 we outline the methods used to 
measure both the velocity divergence power spectrum and 
the matter power spectrum as a function of scale. 

2.1 N-body simulations 

We use the N-body simulations carried out by Li et al. (2012, 
2013). These simulations were performed using a modified 
version of the mesh-based A-body code RAMSES (Teyssier 
2002). Assuming a ACDM cosmology, the following cosmo¬ 
logical parameters were used in the simulations: = 0.24, 

fl DB = 0.76, h = 0.73 and a spectral tilt of ns = 0.961 (in 
agreement with e.g. Sanchez et al. 2009). The linear theory 
rms fluctuation in spheres of radius 8 h _1 Mpc is set to be 
os = 0.769. The simulations use A = 1024 3 dark matter 
particles to represent the matter distribution in a compu¬ 
tational box of comoving length 1500h -1 Mpc. The initial 
conditions were generated at 2 = 49 using the MPgrafic 1 

1 http://www2.iap.fr/users/pichon/mpgrafic.html 
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code. The errors on the power spectra in this work are cal¬ 
culated from the variance in the two point statistics from six 
simulations of the same cosmology initialized with different 
realizations of the dark matter density field. 

We use the publicly available halo catalogues from the 
MultiDark simulation (Riebe et al. 2011; Prada et al. 2012) 
which has a computational box size of L = 1000h _1 Mpc on 
a side. These halos have been identified using the Bound- 
Density-Maxima algorithm (Klypin & Holtzman 1997). The 
halo sample we use in this work consists of all haloes with 
M ^ 5 x 10 12 /i -1 Mq at 2 = 0 and 2 = 1. The error on 
the halo power spectrum in a spherical shell of width 5k is 
estimated using the following formula derived by Feldman 
et al. (1994): 


<7 

P 


( 2^) 2 ( 
Vk 2 5k V 



( 1 ) 


where h is the number density and V is the volume. We 
measure the linear bias, b, for this sample of halos by fit¬ 
ting to the ratio b = -J{5n5p )/{<5lt<5£t) on large scales 
k < 0.1/iMpc -1 , where 5h is the nonlinear halo overden¬ 
sity in Fourier space. Here (<5lt<5lt) is the z = 0 linear the¬ 
ory power spectrum generated using CAMB with the same 
cosmological parameters used in the MultiDark simulations. 


2.2 Measuring the density and velocity fields 

The nonlinear matter and halo power spectra are measured 
from the simulations by assigning the particles to a mesh 
using the cloud in cell (CIC) assignment scheme (Hockney 
& Eastwood 1988) onto a 512 3 grid and performing a fast 
Fourier transform (FFT) of the density field. To compensate 
for the mass assignment scheme we perform an approximate 
de-convolution following Baumgart & Fry (1991). 

Measuring the velocity divergence field accurately from 
numerical simulations on small scales can be difficult if a 
mass weighted approach is used as in Scoccimarro (2004); 
Pueblas & Scoccimarro (2009); Jennings et al. (2011). Some 
volume weighted measures of the velocity field have also 
been developed (see e.g. Bernardeau & van de Weygaert 
1996; Colombi et al. 2007) including the Delaunay tessella¬ 
tion field estimator (DTFE) method (Schaap 2007; Cautun 
& van de Weygaert 2011). 

In the mass weighted approach, simply interpolating 
the velocities to a grid, as suggested by Scoccimarro (2004), 
gives the momentum field which is then Fourier transformed 
and divided by the Fourier transform of the density field, 
which results in a mass weighted velocity field on the grid. 
One of the main problems with this approach is that the 
velocity field is artificially set to zero in regions where there 
are no particles, as the density is zero in these empty cells. 
Pueblas & Scoccimarro (2009) also found that this method 
does not accurately recover the input velocity divergence 
power spectrum on scales k > 0.2/iMpc^ 1 interpolating the 
velocities of 640 3 particles to a 200 3 grid. Using simulations 
of 1024 3 particles in a 1.5ft _1 Gpc box, Jennings et al. (2011) 
found that the maximum grid size that could be used was 
350 3 without reaching the limit of empty cells. 

In this paper the velocity divergence fields are mea¬ 
sured from the N-body simulations using the DTFE method 
(Schaap 2007; Cautun & van de Weygaert 2011). This code 
constructs the Delaunay tessellation from a discrete set of 


points and interpolates the field values onto a user defined 
grid. For the L box = 1500/V 1 Mpc simulation we generate 
all two point statistics on a 512 3 grid. We have verified that 
our results do not change when we increase the grid size to 
1024 3 , demonstrating that our two point clustering statis¬ 
tics have converged on the relevant scales in this paper. The 
velocity divergence field is interpolated onto the grid by ran¬ 
domly sampling the field values at a given number of sample 
points within the Delaunay cells and than taking the aver¬ 
age of those values. The resolution of the mesh used in this 
study means that mass assignment effects are negligible on 
the scales of interest here. Throughout this paper the ve¬ 
locity divergence is normalized to a dimensionless quantity 
6 = —V ■ v/(aH), where v is the peculiar velocity, H is the 
Hubble parameter and a is the scale factor. 

It has recently been shown that there exists a non- 
negligible velocity bias on large scales between the halo 
and dark matter velocity fields. This is a statistical man¬ 
ifestation of sampling effect which increases with decreas¬ 
ing number density (see e.g. Biagetti et al. 2014; Zheng 
et al. 2014; Jennings et al. 2015). We use haloes of mass 
M ^ 5 x 10 12 h~ 1 M@ from the MultiDark simulation which 
have a number density of n = 1.23x 1CU 2 (Mpc//i)~ 3 at 2 = 0 
so that the velocity bias is negligible on the relevant scales 
discussed in this paper. Note that certain methods of mea¬ 
suring either the velocity or velocity divergence field, e.g. the 
nearest grid point method, can induce extra sampling effects 
in addition to the statistical bias mentioned above (see e.g. 
Zhang et al. 2014), the DTFE method does not suffer from 
the same sampling effects (Schaap 2007) and will not impact 
our analysis which is restricted to scales k < 0.45/iMpc _1 . 


3 THE DENSITY - VELOCITY FIELD 
RELATION 


3.1 Linear theory 


At large scales the Universe is homogeneous and the fluctu¬ 
ation fields 5(x, t) = p(x,t)/p — 1, v(x, t), <f>(x, t) are small 
compared to the smooth background contributions. An Eu- 
lerian approach to density fluctuations relies on a truncation 
of the full Vlasov equation and the imposition of an equa¬ 
tion of state. Under the assumption of zero shear, the linear 
regime is then described by the continuity and Euler equa¬ 
tions, 

^l + V.v(*, T) = 0 (2) 

)V ) .:// v;y , i = -V$(x,t), (3) 

where df = adr. The linear theory growth rate, fur is de¬ 
fined as the logarithmic derivative of the overdensity field, 
and is dependent on the cosmological parameters, 


fur(tlm, Ua) 


dlnd 

dlna 


(4) 


The growing mode solution for 5(x, t) admits a product form 
in which it separates as 5{x, t) = D(t)5(x, 0), where D is the 
linear growth factor. For this product form the linear growth 
rate becomes the logarithmic derivative of the growth factor, 
fur(t) = dlnD(f)/dlna. 

Together with the linear continuity equation, we find 
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fc(/iMpc : ) fc(/iMpc x ) 

Figure 1. The upper and lower left panels show the ratio of Re(0(k)/5(k)) (Im(0(k)/<5(k))) as a function of wavenumber, k, measured 
from the simulations at z = 0. The joint pdf 'P(log|0| 2 , log|<5| 2 ) is plotted in the top right panel. The ratio of the magnitudes |0| 2 /|5| 2 
for each mode is plotted in the lower right panel as a function of scale. In all panels the linear theory prediction for the 6 — 8 relation is 
shown as a red dashed line. 



that the velocity divergence and over density fields are simply 3.2 A nonlinear stochastic relation between the 

related as density and velocity fields 


fl(x,t) := V = -/L T (n m ,n A )(5(x,t) (5) 

where we define 9 as the velocity divergence in units of (aH) 
and v(x,t) is the comoving peculiar velocity. Since we are 
within the linear regime, this relation carries over trivially to 
Fourier space, where 6( k, t) = — /i,T<$(k, t). Put another way, 
the linear regime is special in that it admits the introduction 
of a linear growth rate a) that is independent of 

the scale at which we measure the perturbations. 

However, we do not expect this relation to hold once 
the density fluctuations in the fields become large, and non¬ 
linear growth starts to generate mode-coupling. In what fol¬ 
lows we shall analyse to what extent it is possible to sensibly 
extend the central relation Eq. (5) beyond the linear regime, 
and to provide meaningful insights into bulk characteristics 
that arise from non-linearities. We find that the relation is 
modified in essentially two ways: firstly one finds a growth 
factor that is scale dependent due to non-linearities, and 
secondly we find that the deterministic one-to-one relation 
between 6 and 5 is weakened to a stochastic relation. In sec¬ 
tion 4.4 we describe how these results are understood from 
the perspective of perturbation theory, and in section 5 dis¬ 
cuss implications for redshift space distortions. 


Random fields in cosmology are used to represent a single 
realization of the dark matter distribution within a given 
cosmology. As these fields evolve under gravity, nonlineari¬ 
ties give rise to a growth in structure which induce correla¬ 
tions between different scales. The full nonlinear equations 
of motion in Fourier space are given by 

-^jjd T S(k, t) + 9( k, r) as — J d 3 kiA(ki, k - ki)0(ki)<S(k - ki) 

(6) 

d T 0{ k, r) + aH9{ k, r) 

+ ^Q m aH5{ k, t) = — j d 3 kifl(ki, k - ki)0(ki)0(k - ki), 

( 7 ) 


where we have the mode-coupling functions 


A(ki, k 2 ) 


B(ki,k a ) 


(k t + kg) ■ ki 

k\ 

|ki+k 2 | 2 (ki-k 2 ) 

2k\k,2 


( 8 ) 


The terms on the right hand side of both Eq. (6) and 
(7) encode the nonlinear evolution of the fields (see e.g. 
Bernardeau et al. 2002, for a review). Computing the per- 


© 0000 RAS, MNRAS 000, 000-000 























Nonlinear stochastic growth rates and redshift space distortions 5 


turbative components for these nonlinear contributions is 
an on-going challenge, and the complexity increases rapidly 
with higher order terms (see e.g. Crocce & Scoccimarro 2006, 
2008). In what follows we construct a phenomenological ap¬ 
proach to describe the breakdown of the linear theory rela¬ 
tionship in simple terms. 

As mentioned, the departure point from linear theory 
that provides our focus is the relationship between the over¬ 
density field and the velocity divergence. To illustrate this 
breakdown the upper (lower) left panel of Figure 1 shows the 
scatter in the ratio of Re[0(k)/5(k)] and Im[0(k)/<5(k)] as a 
function of wavenumber, k, measured from the simulations 
at redshift z = 0. 

Significant scatter exists in the Fourier modes about the 
linear theory relation (red dashed line), and which increases 
as a function of scale. In the lower right panel of this figure 
we also plot the ratio of the magnitudes |0(k)|/|<5(k)| as a 
function of scale, which demonstrates that this scatter is not 
due to an arbitrary phase differences between the modes, 
and which could have cancelled when computing the two 
point statistics of the fields. The scatter in the 9 — 8 relation 
is also shown in the upper right panel of Fig 1, where we 
plot the PDF of the logarithm of \9 2 \ and |<5 2 |. 

Firstly, we introduce a conditioned velocity divergence 
quantity (6* | <5) that is derived from the conditional dis¬ 
tribution P(#|<5). More explicitly, we define (#(k,f)|<5) := 
f T>8 P(#|(5)#(k, f) for the conditional expectation value of 9. 
The resultant term has a dependence on the particular over¬ 
density that is realised. In the linear regime a direct relation 
exists between 8 and 8, and corresponds to a delta function 
distribution in P(9,8) for which P(#|<5) is perfectly sharp, 
or “deterministic”, and encodes the relation 9 = —f lt8. 
However we can extend this to a more general scenario that 
drops this sharp relation in favor of a stochastic one. We 
define a growth rate /i(n m , Da, k), in momentum space, as 

ft (n m ,H a, k) := -^y<*(k, t)\S) els . (9) 


Here the generalized growth rate now has an explicit de¬ 
pendence on the overdensity field that is being conditioned 
on, in addition to a potential scale dependence. Importantly, 
in the linear regime this function coincides with the linear 
growth rate Jlt, but more generally becomes a stochastic 
quantity for which moments can be computed. 

To estimate the non-linear distortions to the effective 
growth rate, it is instructive to compute the following mo¬ 
ments 

f ... mm —(fs8 2 )s 
J ' ( 8 2 ) ( 8 2 ) y 1 

p. (mm) (fis 2 ) s 

J ■ ( 8 2 ) (< 5 2 ) ’ ^ ’ 


where by definition, (8) = 0. Here (■)$ denotes an ensemble 
average with respect to the probability distribution function 
P(S), however from now on we will omit the subscript from 
any ensemble average notation, for simplicity. In the linear 
regime we automatically have that |/| = |/| = /LT(D m , Da), 
as expected. 

In addition to these non-linear distortions to /, we re¬ 
call that the essential connection between 8 and 8 gradually 
becomes diluted to a stochastic one. This can be quantified 
through the fluctuations of #(k) about the conditional ex¬ 


pectation. In particular, we consider the following random 
held 


a(M):=0(k,i)-(0(k,i)|*> (12) 


whose variance provides a suitable measure, and is given by 


2 = <^> = (8 2 ) - (m 2 ) 

“ (<5 2 ) ( 8 2 ) 


(13) 


Again if the linear continuity equation, Eq. (5), holds then 
a a = 0 and stochastic relation vanishes. Section 4 contains 
a closer examination of the nonlinear moments in Eqs. (10) 
and (11), and the magnitude of a measured from the N-body 
simulations as a function of both scale and redshift. 

More generally it is seen that the quantity / is related 
to the expected velocity divergence at a particular scale 
through the relation 

(9) = — J d 3 fciA(ki, k — ki)x (14) 

[/(J(kr)J(k - kr)> + (a(kr)«5(k - kr)>] 


which follows from the full non-linear continuity equation. 
A parallel relation for / can be obtained, and from the Euler 
equation we End 


(8) = -Jd 3 ki (15) 

[P(ki, k - k 1 )(/ 2 (,f(k 1 )5(k - kO) + (a(ki)a(k - ki)» 

—A(ki, k - ki)(/(5(ki)5(k - ki)) + (a(ki)i(k - ki)»] , 

where a quantifies the deviation from a deterministic re¬ 
lation between 9 and 8 and we have used the fact that 
(0(k 2 )0(k;O) = «6>(k 1 )|5)(0(k 2 )|<5» - <a(k 1 )a(k 3 )>. 

The quantities f and / are readily extracted from sim¬ 
ulations, for which we restrict the analysis of velocity and 
overdensity fields to large scales to avoid issues associated 
with the measurement of the velocity held in an unbiased 
way (Pueblas & Scoccimarro 2009; Jennings et al. 2011). 

As a side note, it is common to address the de¬ 
gree of stochasticity between two random functions X 
and Y by measuring the cross correlation coefficient r = 
{XYi/y/wmm as a function of scale. This is a dif¬ 
ferent notion of stochasticity to the one discussed in this 
paper, and relates to either a bias between the two helds 
at the level of the two point functions, (|T| 2 > = &2 Pt (m 2 ) 
or a more specific local bias Y = bi OC a.\X. As pointed out 
by Dekel & Lahav (1999), the bias between the two point 
statistics follows from a local deterministic bias, and is the 
square of the local bias but the converse does not necessar¬ 
ily follow. In this case the cross correlation coefficient is a 
measure of r = &2 P t/&iocai and is not necessarily unity. Here 
the bias &i oca i could represent the familiar bias between the 
mass and halo/galaxy overdensity or we could view it as the 
linear growth rate in the overdensity - velocity divergence 
relation in linear theory. As pointed out by Seljak & Warren 
(2004) the cross correlation coefficient can be close to unity 
despite huctuations about a local bias being large. 


3.3 Decomposition of two-point functions 

It is also instructive to decompose the two-point func¬ 
tions (9(ki)8(k2)) and ( 6(k\)9{k2 )) into contributions com- 
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fc(hMpc *) 

Figure 2. The first moments of / and /, given in Eqs. (10) and 
(11), measured from the simulations at z = 0 are shown as a 
blue dot dashed and black dashed line respectively. The linear 
theory growth rate is shown as a red solid line in this figure. 
The variance in / measured from six realizations of the same 
cosmology is shown as the grey shaded region. 

ing from the nonlinear corrections and stochasticity in the 
0—8 relation. The two-point function for a(k, t) in Eq. (12), 
decomposes as 

(a(k,t)a*(k',t)) = ((0(k,t)|5>(r(k , ,t)|,5»-(0(k,f)r(k',i)}. 

From this the two-point functions of interest - the auto and 
cross power spectra between the conditional mean of (#|<5) 
and 5 - can be expressed as 


(9i 5* 2 ) = h 2 (5i5* 2 ) + (ax5* 2 ) 

(16) 

= <(0r| + 


<M;> = /? 2 <*i£) + <aia3> 

(17) 


= «0i|.S)(02|a» + (ona 2 ). 

Here we employ the short-hand notation Xi for one-point 
quantities -Y(ki) and Y iy for two-point quantities V'(k,;. k,). 
Note that / and / are now evaluated as two-points functions. 
This separates out the nonlinear and stochastic components, 
as defined in Section 3.2 in a natural way, and empha¬ 
sizes the different dependence on stochasticity for the auto¬ 
correlation and cross-correlation spectra. Also note that this 
approach is in contrast to previous studies (Scoccimarro 
2004; Percival & White 2009; Jennings et al. 2011; Jen¬ 
nings 2012; Crocce et al. 2012; Jennings et al. 2012) which 
compare the ensemble averaged statistics Pgs = (9 i5 2 ) and 
Pgg = ( 6 x 62 ) with Pss = (Ji $ 2 ) as a function of scale. In Sec¬ 
tion 4 we present the measurements of these two point func¬ 
tions and test the decomposition into nonlinear and stochas¬ 
tic components given in Eqs. (16) and (17). 

4 RESULTS 

We now provide a more detailed account of how the quan¬ 
tities introduced in the previous section behave in practice. 


The values of / and / are computed in Section 4.1 at differ¬ 
ent scales, and compared with linear theory. The deviation 
of 9 from the conditional mean (#|<5) is addressed in Section 

4.2, both as a function of scale and redshift. We verify that 
the decomposition of the two point statistics into nonlinear 
and stochastic parts as defined in Section 3.2 is reproduced 
within the simulation and we present the measured two point 
statistics in each case. 

It turns out that halos display these features more dra¬ 
matically than dark matter, and this is discussed in Section 

4.3, where we measure (#|d) and the associated two point 
functions for halos with masses M ^ 5 x 10 12 Mgft _1 from 
the MultiDark simulations. An obvious question is: to what 
degree are these features reproduced by existing perturba¬ 
tive results, and do the decompositions presented simply cor¬ 
respond with a particular perturbative order? To this end, 
in Section 4.4 we compare our results with standard pertur¬ 
bation theory to third order and second order Lagrangian 
Perturbation Theory (2LPT) predictions for the two point 
functions (69) and (95). 

4.1 Nonlinear growth functions / and / 

The degree to which the moments given in Eqs. (10) and 
(11) in Section 3.2 differ from /lt are a measure of the 
deviations from linearity, and provide effective non-linear 
growth rates. In Fig. 2 we plot these two moments, / and /, 
as a red solid (for /lt), blue dot dashed and black dashed 
lines respectively, measured from the non-linear dark matter 
density field in the simulations at 2 = 0. 

Note that the two moments / and / that are plotted are 
the average of six N-body simulations initialized with dif¬ 
ferent realizations of the matter density field at early times. 
The variance of / from these six simulations is shown as a 
grey shaded region. 

We find a notable difference between the three growth 
rates, and even on large scales, such as k < O.lhMpc -1 , nei¬ 
ther / nor / correspond to the linear theory growth rate /lt- 
We find that the ratio of the two-point functions y/ (99) / (55) 
and (59) / (55) do converge to the linear theory result /lt on 
much larger scales k < 0.02/iMpc _1 in agreement with pre¬ 
vious work (Scoccimarro 2004; Percival & White 2009; Jen¬ 
nings et al. 2011; Jennings 2012; Crocce et al. 2012; Jennings 
et al. 2012). Taking the decomposition of each of these two 
point functions given in Eqs. (16) and (17) into account this 
implies that on large scales the ratio of the stochastic two 
point functions ( a5 ) and (aa) to (55) is at a minimum 10- 
15% of /lt at k < O.lhMpc -1 . We demonstrate that both 
of these decompositions are valid in Section 4.2. 

4.2 The stochastic relation between 9 and 5. 

In Fig 3 we plot the conditional expectation (#|<5) as a or¬ 
ange dashed line. This is the average over six realizations, 
measured from the simulations by simply binning in Red(k) 
and finding the mean Re#(k) at 2 = 0 at the two scales 
k = O.lhMpc^ 1 (left panel) and k = 0.2/iMpc _1 (right 
panel), while the linear theory relation between 5 and 9 is 
plotted as a red line. The real component of the Fourier 
modes measured from one simulation at each wavenumber 
are shown as black dots. At each scale k there is signifi¬ 
cant scatter between the Fourier modes measured from the 
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Figure 3. The conditional expectation (0(k)|(5) (orange dashed line) at k = 0.1/iMpc -1 (left) and k = 0.2hMpc -1 (right) at z = 0, 
together with the linear theory relation between 5 and 6 (red line) and the real Fourier modes (black dots) measured from a ACDM 
simulation (/-box = 1500Mpc/h) at z = 0. In the right panel we show the conditional expectation (0|<5)lpt from second order Lagrangian 
perturbation theory as a cyan dot dashed line. The inset panel shows the ratio of the two point function ((0|<5)(0|<5))/{00)lpt measured 
from the simulations at z = 0 as a function of scale as a blue solid line. 



300000 


200000 - 


100000 


0 - 


-100000 


-200000 - 


T-T 

k= 0.2hMpc _1 

z - 0.4 



-300000 


- linear theory 

- ■ (0|«5> 

J_I_L 


_L 


-200000 -100000 0 

Re Ml 


_L 


100000 200000 


Figure 4. The conditional expectation (0(k)|<5) (orange dashed line) at k = 0.1/iMpc -1 (left) and k = 0.2/iMpc —1 (right) at z = 0.4, 
together with the linear theory relation between 5 and 6 (red line) and the real Fourier modes (black dots) measured from a ACDM 
simulation. 


simulations, and (6* | <5) differs from linear perturbation the¬ 
ory predictions of 6 = —J'ltS. It is also clear from these two 
panels that the difference between (#|5) and — furS increases 
with increasing k into the nonlinear regime. 

There are two notable effects which are evident from Fig 
3. Firstly, the nonlinearity we are describing with the con¬ 
ditional mean {#|<5} manifests as an approximate rotation 
about the linear theory prediction (orange dashed line in 
Fig. 3 compared to the solid red line) which is linear in S but 
with a scale dependent coefficient i.e. (0|5) ~ —fuTS + c(k)5, 


where c is an increasing function of scale. In Section 4.4 we 
show that this functional dependence can be explained on 
large scales by second order Lagrangian perturbation the¬ 
ory. The second thing to note from these plots is that the 
stochastic scatter around (#|<5) is nonzero and increases with 
increasing wavenumber k. At both scales we find that for 
5 > 0(< 0) the mean relation (#|<5) is larger (smaller) then 
the linear theory prediction, corresponding to an effective 
growth factor that is larger than linear theory. 

The corresponding plot at 2 = 0.4 is shown in Fig. 4 
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Figure 5. The ratio of the two point functions in Eqs. (16) 
and (17) to ( 61 S 2 ) and [0\ 82 ) measured from the simulations at 
z = 0 are shown as red dashed and green dot dashed lines respec¬ 
tively. The ratios of the two point functions ([0\ 1 ^) 62 )/ (0\ 62 ) and 
{(01 |(5) (02 |<5))/(0i 62 ) are shown as a blue dot and black dashed 
line respectively. The shaded cyan and grey regions show the vari¬ 
ance of these ratios measured from six simulations. 
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for the same two scales. At higher redshifts we see the same 
trend with (0|<5) behaving as a rotation away from — /lt< 5 
for linear theory. This difference increases with increasing 
wavenumber although this difference is smaller then at z = 
0 due to increased nonlinear growth at later redshifts as 
expected. We also note a reduction is the range of 5 values 
at z = 0.4 compared to z = 0. An identical number of modes 
have been used at each scale and redshift. 

The decomposition of the two point functions ( 95) and 
( 86 ) into nonlinear and stochastic parts, as in Eqs. (16) and 
(17), is readily verified numerically. In Fig. 5 we plot the 
ratios of the RHS of Eqs. (16) and (17) to ( 8162 ) and ( 81 S 2 ) 
as red dashed and green dot dashed lines respectively. We 
find these ratios are unity which verifies the decompositions 
in Eqs. (16) and (17) from the simulations. This result is 
non-trivial as all of the two point statistics have been mea¬ 
sured independently from the simulations i.e. ((0|<5)<S) is an 
ensemble average over the mean 8 given 5 (orange dashed 
line in Fig. 3) times 5. This is in contrast to (85) which is 
the ensemble average over each 8 and 5 (black dots in Fig. 

3). 

In Fig 5 we also plot the ratios ((8\5){8\5))/( 88 ) and 
((0|<5)<5)/(0<5) measured from the simulations at z = 0 as 
black dashed and blue dot dashed lines respectively. The 
shaded regions in this plot represent the variance amongst 
six realizations of the same cosmology. We find that the 
stochastic components contribute ~ 10% to the two point 
function ( 8182 ) while it contributes about a 1% to ( 8162 ) at 
k < 0.2fcMpc -1 . The stochastic component of the velocity 
divergence auto power increases to approximately 25% by 
k = 0.45hMpc -1 . In the upper and lower panels of Fig. 6 
we show similar power spectra ratios to those in Fig. 5 at 
z = 0.4 and 2 = 1 respectively. We find that the stochas- 
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Figure 6. The ratio of the two point functions in Eqs. (16) and 
(17) to ( 0 i< 52 ) and (0\ 62 ) measured at at 2 = 0.4 (top panel) 
and at z = 1 (lower panel) are shown as red dashed and green 
dot dashed lines respectively. The ratios of the two point func¬ 
tions (( 0 i 1 5) 62 )/ (0i (52) and ((011<5) (021<5))/(0102) are shown as a 
blue dot and black dashed line respectively. The shaded cyan and 
grey regions show the variance of these ratios measured from six 
simulations. 
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tic component of the velocity divergence power spectrum is 
slightly reduced at higher redshifts as there is less nonlinear 
growth present at earlier times which would induce a larger 
variation in 8 from (0|<5). 


4.3 Behaviour of (0|<5) and a for dark matter halos 

In Fig. 7 we show the conditional expectation (0(k)|<5) as 
a purple dashed line, measured at k = O.lhMpc - (left 
panel) and k = 0.2/iMpc^ 1 (right panel) from the MultiDark 
simulations using halos with masses M < 5 x 10 12 . 

The real Fourier modes R ednik) and ReJsr(fc) are shown 
as grey dots in both panels. We plot the linear theory pre¬ 
diction 8 h = —fur/bSH, where b is the linear bias on large 
scales as a red solid line. The red dotted lines either side 
of the linear theory prediction represent a ±10% error in 
the linear bias. For this halo sample we find that the linear 
bias is b ~ 0.81 ± 0.09 and is reasonably linear on scales 
k ^ 0.2/iMpc _1 (see also e.g. Jennings et al. 2015). 

It is clear that there is significant scatter about the 
mean (0(k)|<5) and that this conditional expectation differs 
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Figure 7. The conditional expectation (0(k)|<5) (purple dashed line) measured using halos with M < 5 X 10 1 21 in the MultiDark 
simulation at k = O.lhMpc” 1 (left) and k = 0.2 /AT pc 1 (right) at z = 0. The linear theory relation Oh = — f\;v /hi) h , where b is the 
linear bias on large scales is shown as a red solid line. The red dotted lines either side of the linear theory prediction represent a ±10% 
error in the linear bias. The real Fourier modes measured using this halo catalogue at each wavenumber are shown as grey dots. 
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Figure 8. The ratio of the two point function 
(( 011 <5)(<?2|<5))/(0 i02 ) measured at z = 0 and 2 = 1 using 
halos with M < 5 X 10 12 Mg/i -1 from the MultiDark simulation 
are shown as black dashed and blue dot dashed lines respectively. 
The shaded tan region represents the error on the measured 
power spectra given in Eq. (1) in Section 2.1. 
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ter (Multidark) simulations, which may account for some of 
these differences. 

In Fig. 8 we plot the ratio of the two point function 
{{ 6 i\ 5 ){ 62 \ 5 ))/{ 6162 ) measured at 2 = 0 and 2 = 1 for the 
same halo sample as black dashed and purple dot dashed 
lines respectively. The shaded tan region represents the er¬ 
ror on the measured power spectra given in Eq. (1) in Section 
2.1. We have also verified that the decomposition of the two 
point functions into nonlinear and stochastic parts, as given 
in Eqs. (16) and (17), holds for the halo two point functions. 
We have omitted this from Fig. 8 for clarity. Therefore any 
deviation from unity in this figure indicates the magnitude of 
the stochastic component. We find that the stochastic com¬ 
ponent of the two point function ( 6162 ) is significant and ap¬ 
proximately a constant fraction (~15%) at k < O^SftMpc” 1 
at both 2 = 0 and 2 = 1. The differences between the halo 
sample and the dark matter distribution, in how the condi¬ 
tional mean deviates from the linear theory predictions and 
the scatter around that mean as a function of wavenumber, 
cannot be only due to a difference in cosmological parame¬ 
ters. As shown in Fig. 8 we find that the velocity divergence 
auto power spectrum has a larger stochastic component in 
the halo sample compared to the dark matter at both red- 
shifts. 


from the linear theory prediction by an approximate rota¬ 
tion. If we compare these results with Fig. 3 in Section 4.2 
for the dark matter we see that at the same redshift, the 
deviation of (0(k)|<5) from the linear theory prediction and 
the scatter about the conditional mean given by a, is much 
larger for the halo sample then for the dark matter at both 
k scales. Note these two simulations have slightly different 
cosmologies, for example = 0.24 (0.27) in the dark mat- 


4.4 The relation to standard and Lagrangian 
Perturbation Theory 

In this section we connect the formalism presented in Sec¬ 
tion 3.2 to perturbation theory methods. First we consider 
standard perturbation theory predictions for both the veloc¬ 
ity divergence auto and cross power spectra, (99) and (95), 
computed up to third order (see e.g. Bernardeau et al. 2002, 
for a review). The nonlinear velocity divergence auto P( k) 
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Figure 9. The ratio of the one loop perturbation theory pre¬ 
diction for the velocity divergence cross and auto power spectra 
(Eq. (18)) to (OS) and (06) measured from the dark matter only 
simulations are shown a grey dot dashed and black dashed line 
respectively. The pink and cyan error bars show the variance in 
these ratios from six simulations with different realisations of the 
initial density field. Similar ratios of Pn + P 22 calculated from 
one loop perturbation theory for the both cross and auto power 
spectrum are shown as blue solid and orange dot dashed lines as 
given in the legend. 


computed from third order perturbation theory is given by 
(flfl)pT(k) = P(k) + P 22 ( k) + 2P i 3 (k) (18) 


where P(k) denotes the linear power spectrum and the scale 
dependent functions P 22 and Pi 3 are given by 


P 22 (k) = 6 P(k) j d 3 gG 3 (k,q)P(q) (19) 

Pl 3 (k) = J d 3 q[G 2 (k - q, q)] 2 P(k - q)P(q), ( 20 ) 

where the kernel G 2 is given by 

c 2 (k .«0 = y+ 5*=-$(- + £) 

+ 1 (kikj - 1 Sij ^ (qiqj - l<5i^ > ( 21 ) 


where fj ,2 = 26/21 and the angle averaged G 3 kernel is given 
in e.g. Eq. (69) in (Scoccimarro 2004). A similar expression 
for the velocity divergence cross power spectrum to third 
order can also be found in (Scoccimarro 2004). 

In order to compare the formalism in this paper, which 
decomposes the velocity divergence two point statistics into 
nonlinear and stochastic elements as given in Section 3.2, 
with perturbation theory methods we simply calculate the 
individual power spectra in Eq. (18) and compare them with 
the measure two point velocity divergence statistics. In Fig 
9 we show the ratio of the one loop perturbation theory 
predictions for the velocity divergence cross and auto power 
spectra to ( 68 ) and (96), measured from the dark matter 
only simulations at z = 0 , as a grey dot dashed and black 
dashed line respectively. The pink and cyan error bars show 


the variance in these ratios from six simulations with dif¬ 
ferent realizations of the initial density field. The ratios of 
Pn +P 22 calculated from 1 loop perturbation theory for the 
both cross and auto power spectrum to ( 68 ) and ( 66 ) are 
shown as blue solid and orange dot dashed line respectively. 

By comparing Figs. 5 and 9 we can see that at the level 
of 3rd order perturbation theory that the standard pertur¬ 
bation theory prediction and the formalism in this paper 
deviate substantially and no simple identification can be 
made. Even on large scales, k < 0.05/iMpc^ 1 , where the 
perturbation theory predictions match the measured power 
spectra from the simulations to ~ 5% we cannot simply re¬ 
late the mode coupling terms Pi 3 , which are negative, to the 
stochastic power spectra (aa) and (a<5). 

Next we consider second order Lagrangian perturbation 
theory (2LPT) predictions for the 6 — 8 relation (Gramann 
1993; Bouchet et al. 1995; Bouchet 1996; Melott et al. 1995; 
Kitaura et al. 2012). Lagrangian perturbation theory rep¬ 
resents a alternative framework to the Eulerian approach, 
and the non-linear analysis is based around the trajectories 
of individual fluid elements. Of central importance is the 
displacement field 'F(q), which provides a mapping from 
initial Lagrangian coordinates q to final Eulerian coordi¬ 
nates x given by x(r) = q+ 4t(q, r) (see e.g. Bouchet 1996; 
Bernardeau et al. 2002, for a review). The linear solution 
for the equations of motion coincide with the Zel’dovich ap¬ 
proximation, 

V q -\t (1) = -P(r)5 ( 1 ) (q), (22) 

where 5^(q) is the linear density field and D is the linear 
growth factor normalized to unity at z = 0. In contrast, the 
second order correction to the displacement field (see e.g. 
Melott et al. 1995) takes into account tidal gravitational 
effects as 

v q ■ 4> (2) = \d 2 Y, (*$»& - [*Sf) , (23) 

where = d^i/dqj and D 2 is the second order growth 
factor given by D2 ~ —3/7 D 2 fim 143 . The Lagrangian poten¬ 
tials cjJ 1 ^ and cf>^ are defined such that V 2 0 ^(q) = 5^(q) 
and V 2 cp( 2 \ q) = <^ 2 ) (q). The 2LPT expressions for the po¬ 
sition become 

*(q) = q PV,^ (1) + P 2 V,</ (2) , (24) 

while the dimensionless velocity divergence is given by 

6 = -P/ lt V 2 /. (1) + P 2 / 2 V 2 0 (2) , (25) 

where /lt is the linear theory growth rate and / 2 = 
dlnP 2 /din a is the logarithmic derivative of the second-order 
growth factor, / 2 « 2 Sim 11 . The gradient terms are given by 

V 2 0 (1) =5 ( 1 ) (q) (26) 

v^ {2) = e(^M-[^] 2 ). ( 27 ) 

i>j 

where (jJJ = d 2 (j> / dqidqj. These equations imply that given 
a robust estimate of the linear overdensity then we can 
obtain a corresponding non-linear velocity divergence 6 from 
2LPT. 

To compare the conditional mean (0|5) measured from 
the simulations in Fourier space with (#|5)lpt, where 6 is 
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the Fourier transform of the 2LPT prediction in Eq. (25), we 
estimate the linear density field 5^ from the nonlinear mat¬ 
ter field 5 m in real space as given by Neyrinck et al. (2009), 
<5 (1) = log(l + 5m) — (log(l + 5m))- The quantity / 2) (k) 
can be obtained by Fourier transforming (k),ij into real 
space, computing the sum and then transforming back. Al¬ 
ternatively one can Fourier transform Eq. (27) directly, and 
then obtain the total 2LPT prediction for #(k). 

In the right panel of Fig. 3 we plot the conditional ex¬ 
pectation (#|<5)lpt, which has been evaluated by the same 
method as described in Section 4.2 using the 2LPT predic¬ 
tion for #(k), as a cyan dot dashed line. From this figure it 
seems that the rotation of the conditional mean away from 
the linear theory prediction is captured well by 2LPT at 
k = 0.2hMpc -1 . In the inset panel of Fig. 3 we also show 
the ratio of the two-point function ((#|5)(#|<I))/(##)lpt mea¬ 
sured from the simulations at 2 = 0 as a function of scale as 
a blue solid line. On scales k < 0.1/iMpc -1 we can see that 
the ratio is very close to unity indicating that the nonlin¬ 
ear effects in the formalism of this paper can be described 
by 2LPT which incorporate the effects of tidal gravitational 
fields on large scales. Recall that the full two point func¬ 
tion {99) can be written as sum of nonlinear and stochastic 
components. Our results indicate that the stochastic com¬ 
ponent, which is approximately 15% of the velocity diver¬ 
gence auto power spectrum at z = 0 and k = 0.1/iMpc ^ 1 is 
not described by the predictions of 2LPT. At smaller scales 
k > 0.1/iMpc -1 we find that the 9-5 relation is still well- 
described by the combined action of a non-linear rotation 
together with stochastic spread, however the predictions of 
2LPT no longer adequately describe the regime. 


5 REDSHIFT SPACE DISTORTIONS 

We begin in Section 5.1 by briefly reviewing the theory of 
redshift space distortions and models that depend on the 
linear growth rate /lt and which are currently in use. In 
Section 5.2 we highlight the problems associated with having 
a well defined notion of the linear growth rate in a redshift 
space distortion model in the presence of a nonlinear and 
stochastic 9 — 5 relation. 

5.1 Redshift space distortion models 

Inhomogeneous structure in the Universe induces peculiar 
motions which distort the clustering pattern measured in 
redshift space on all scales. This effect must be taken into 
account when analyzing three dimensional datasets that use 
redshift to estimate the radial coordinate. Redshift space ef¬ 
fects alter the appearance of the clustering of matter, and 
together with nonlinear evolution and bias, lead the mea¬ 
sured power spectrum to depart from the simple predictions 
of linear perturbation theory. The comoving distance to a 
galaxy, s, differs from its true distance, x, due to its pecu¬ 
liar velocity, v(x) (i.e. an additional velocity to the Hubble 
flow). The mapping from redshift space to real space is given 
by 

s = x + u z z, (28) 

where u z = v ■ z/(aH) and H(a) is the Hubble parameter. 
This assumes that the distortions take place along the line 


of sight, denoted by z, and is commonly referred to as the 
plane parallel approximation. 

On small scales, randomised velocities associated with 
the motion of galaxies inside virialised structures reduce the 
power. The dense central regions of galaxy clusters appear 
elongated along the line of sight in redshift space, which 
produces the ‘fingers of God’ effect in redshift survey plots. 
For growing perturbations on large scales, the overall ef¬ 
fect of redshift space distortions is to enhance the clustering 
amplitude. Any difference in the velocity field due to mass 
flowing from underdense regions to high density regions will 
alter the volume element, causing an enhancement of the 
apparent density contrast in redshift space, 5 s (k), compared 
to that in real space, 5 r (k) (see Hamilton 1998, for a review 
of redshift space distortions). 

Assuming the line of sight component of the peculiar 
velocity is along the 2 -axis, the power spectrum in redshift 
space is given by (Scoccimarro 2004) 

5o(k) + Ps(k) = f'-^e~ itp {e ik * v {l + <5 g (5?)][l + «„(/)]> 

(29) 

where 5 g = b5 is the galaxy overdensity which we shall as¬ 
sume is related by a linear bias, b to the matter overdensity, 
V = u z {x) — u z (x') and f = x — x'. We are also assuming 
that there is no velocity bias between the dark matter and 
galaxies for simplicity. 

Decomposing the vector field into curl and divergence 
free parts, and assuming an irrotational velocity field, we 
can re-write k z u z = ~{k 2 /k 2 )9{k) = —fi 2 9(k) where 9{k) is 
the Fourier transform of the velocity divergence defined in 
Eq. (5). Expanding the exponential term and only keeping 
terms up to second order in the variables 5 and 9, the power 
spectrum in redshift space P s becomes 

5 D {k - k’)Ps(k) = b 2 {5{k)5 *(/?)} - 2 f i 2 b{9{k)5*(k')) 

+ /(0(fc)6>* (£')}• (30) 

If we assume the linear continuity equation holds we can 
re-write this as 

5 D {k-k')P s {k) = {5{k)5*{k'))[b 2 - 26/ LT / + /lt/] 

= S D (k - k')P{k)[b 2 - 2&/ L t/ + /lt/] 

(31) 

which is the Kaiser (1987) formula for the power spectrum 
in redshift space in terms of the linear growth rate /lt, the 
linear bias b and the power spectrum P(k). 

Commonly used models for the redshift space power 
spectrum extend the Kaiser formula by assuming that the 
velocity and density fields are uncorrelated and that the joint 
probability distribution factorizes as P{5,9) = P{9)V{5). 
Examples include multiplying Eq. (31) by a factor which 
attempts to take into account small scale effects, invoking 
either a Gaussian or exponential distribution of peculiar ve¬ 
locities. A popular phenomenological example of this which 
incorporates the damping effect of velocity dispersion on 
small scales is the so-called ‘dispersion model’ (Peacock & 
Dodds 1994), 

P s (k,fj,) = P 3 (fc)(l + fi/) 2 (l — k 2^2 a 2 /2 ) ’ ( 32 ) 

where P g is the galaxy power spectrum, a p is the pairwise 
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velocity dispersion along the line of sight, which is treated 
as a parameter to be fitted to the data and f) = fur/b. 

The linear model for the redshift space power spectrum 
can be extended by keeping the nonlinear velocity power 
spectra terms in Eq. (30). For example Scoccimarro (2004) 
proposed the following model for the redshift space power 
spectrum in terms of Pss, the nonlinear matter power spec¬ 
trum, 

P\k,J) = (Pss(k)+2fi 2 P se (k)+ f i 4 Pee(k)) x e~ (k ^ v)2 , 

(33) 


where Pee = ( 06 ), Pse = (SO) and a v is the ID linear velocity 
dispersion given by 


2 


<7 


V 


1 

3 


Pee(k) 

k 2 


d 3 k. 


(34) 


In linear theory, Pge and Pse take the same form as Pss and 
depart from this at different scales. Using a simulation with 
512 3 particles in a box of length 479/i _1 Mpc, Scoccimarro 
(2004) showed that this simple ansatz for P a (k,fi) was an 
improvement over the Kaiser formula when comparing to 
the results of IV-body simulations in a ACDM cosmology. 
Clearly the inclusion of these nonlinear velocity divergence 
terms gives rise to an improved model of redshift space dis¬ 
tortions in the nonlinear regime. 

In nonlinear models for the power spectrum in redshift 
space there is a degeneracy between the nonlinear bias, the 
difference between the clustering of dark matter and halos 
or galaxies, and the scale dependent damping due to velocity 
distortions on small scales. This degeneracy will complicate 
any measurement of the growth rate using redshift space 
clustering information on small scales. In this work we have 
restricted our analysis of the 9 — 8 relation to large scales for 
the halo sample where the approximation of a linear bias is 
valid. Note also that nonlinearities in the bias between the 
halos and dark matter field affect the J 2 component but not 
the /r 4 coefficient if there is no velocity bias present. 


5.2 Modeling redshift space distortions with a 
nonlinear stochastic 9 — 8 relation 

Firstly, the expansion in Eq. (30) does not assume that 9 
and 8 are uncorrelated (V(S,9) = V(9)V(8)) but instead 
only retains terms which are second order in 9 and 8 . We can 
rewrite Eq. (30) in terms of the main formalism in this paper 
which describes a nonlinear, stochastic relation between 9 
and 8 . Using Eqs. (16) and (17) with the adapted notation 
( 9(k)9(k ')) = (9\92) etc. we can write 

S D {k — k')P s (k) =b 2 (8(k)8*(k')) (35) 

-26 At 2 [((0i|5>5 2 ) + (ai5 2 )] 

+ //[(aia 2 ) + ((0i|<S)(0 2 |<S))]. 

There are a small number of papers that have used per¬ 
turbation theory to find an analytic formula for the condi¬ 
tional mean (#|<5) (see e.g. Chodorowski 1998; Bernardeau 
et al. 1999). Guided by the results in Section 3.2 where (#|<5) 
appears as a rotation from the linear perturbation theory 
prediction which increased with increasing wavenumber k, 
we consider the following simple expression for (#(k)|5) = 
—/LT<f(k) + c(k)<5(k). Putting this into Eq. (35) we obtain 


the expression 

S D (k-k')P s (k) =( 8182 ) (b 2 - 26/r 2 [—/lt + c(k) + 

+/ r 4 [(/ LT -c(k)) 2 + ^i^]). (36) 

A key point that this highlights is that the coefficients 
in front of the J 2 and fi 4 terms no longer have a sim¬ 
ple relation. The receive different contributions from non¬ 
linearity and stochasticity, and cannot be simply written as 
/nlM 2 + /nlM 4 - If the relation between 9 and 8 is determin¬ 
istic ((ai<5 2 ) = 0 and (ai« 2 ) = 0) then, as shown in Section 
4.4, second order Lagrangian perturbation theory provides 
a good description of the nonlinear rotation of the condi¬ 
tional mean (#|<I) away from the linear perturbation theory 
predictions at k < O.lhMpc . The stochastic components 
(ai 82 } and (qiq 2 ) are nonzero at z = 0,0.4 and z = 0 on 
large scales, as can be seen from Figs. 5 and 6, and comprise 
approximately 10% of the velocity divergence auto power 
spectrum on large scales k < 0.1/iMpc -1 . 

It is common practice to try to extract a measurement of 
the linear growth rate, /nl, using the J 2 and /r 4 dependence 
of the measured galaxy power spectrum in redshift space, 
and either the model in Eq. (32) or models which include the 
velocity divergence auto and cross power spectra. If however 
there is a nonlinear and stochastic relation between 9 and 8 
then the correspondence between the coefficients of J 2 and 
/r 4 , and /nl becomes more complex. 

Ideally a perturbative expansion which captures all the 
nonlinearities in Eqs. (6) and (7) would give an accurate 
prediction for the velocity divergence two and higher point 
statistics and their correlations with the matter overdensity. 
This would include the stochastic terms in the formalism in 
this paper which are produced by nonlinear effects. With¬ 
out this exact expansion, it is not straightforward to make 
an explicit connection between the quantity of interest, /nl, 
and parameters in current phenomenological models for two 
point clustering statistics in redshift space, which either as¬ 
sume that 9 and 8 are related by a linear, deterministic re¬ 
lation, or are based on perturbation theory expressions to a 
given order for the 9 — 8 relation. 


6 SUMMARY & CONCLUSIONS 

Up and coming galaxy redshift surveys aim to measure the 
linear growth rate to an accuracy of ~ 1%. This growth 
rate is commonly obtained from a deterministic relation be¬ 
tween the velocity divergence and the matter overdensity 
fields that follows from linear theory. Here we have explored 
a formalism that defines both a nonlinear and a stochas¬ 
tic relation between the velocity divergence and overdensity 
field, 9 = V ■ v(x.,t)/aH and 8 , which is based on an exten¬ 
sion of linear theory to a relation in terms of the conditional 
mean (#|5), together with fluctuations of 9 around this non¬ 
linear relation. 

Using N-body simulations of dark matter particles that 
follow the gravitational collapse of structure over time, we 
measure both the nonlinear and stochastic components and 
verify that this decomposition of the two point clustering 
statistics is reproduced within the simulation. We find that 
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the net effect of the nonlinearity manifests itself as an ap¬ 
proximate rotation of (#|d) away from the linear theory pre¬ 
diction —fur 5, and which increases as a function of scale. 
The scatter about this mean value corresponds to stochas- 
ticity, or variance, of 9 around (#|<5) and which is nonzero 
on all scales. The stochastic contribution to the velocity 
divergence auto power spectrum is approximately 10 % at 
k < 0.2hMpc _1 at z = 0. The stochastic component of the 
velocity divergence auto power increases to approximately 
25% by k = 0.45hMpc _1 . 

We examine two scales in detail, k = 0.1/iMpc -1 and 
k = 0.2/iMpc - , and find that the scatter around the 
mean value ( 0 |<S) is nonzero and increases with increasing 
wavenumber k. At both scales we find that for 5 > 0(< 0) the 
mean relation (#|d) is larger (smaller) then the linear theory 
prediction. We find that both of these trends for the stochas¬ 
tic relation and nonlinearity are visible at higher redshifts, 
s = 0.4 and z = 1 but with a reduced level of stochasticity 
overall due to less nonlinear growth at high redshifts. Using 
a halo sample with M ^ 5 x 10 12 M@/j -1 from the MultiDark 
simulation we find that both the stochasticity and nonlin¬ 
earity in the 6 — 5 relation are larger for halos compared to 
the dark matter. We find that the stochastic component of 
the two point function (99) is significant and approximately 
a constant fraction (15%) at k < 0.25/iMpW 1 at both z = 0 
and z = 1 . 

The relation with perturbative methods was also ex¬ 
plored, and a computation of the velocity divergence auto 
(99) and cross (95) power spectra using one loop stan¬ 
dard perturbation theory reveal that at this level the stan¬ 
dard perturbation theory prediction and the formalism in 
this paper are not equivalent. Even on large scales, k < 
0.05/iMpc -1 , where the perturbation theory predictions 
match the measured power spectra from the simulations to 
~ 5% there is no simple correspondence between the mode 
coupling terms P 13 , which are negative, to the stochastic 
power spectra (a<5) and (aa). 

Using an expression for 9 computed from second order 
Lagrangian perturbation theory we find that the rotation 
of the conditional mean (#|<5) away from the linear theory 
prediction is well described by the conditional expectation 
(#|<5 )lpt from 2LPT on scales k < 0.1/iMpc -1 . This indi¬ 
cates that the nonlinear components in the formalism can 
be described through the inclusion of tidal effects of the 
gravitational field at second order. 

The central features discussed also have an impact on 
the extraction of the linear theory growth rate from models 
of two point functions in redshift space given the level of non¬ 
zero stochasticity which we have measured. It is common 
practice to try to extract a measurement of the linear growth 
rate using the ^t 2 and /x 4 dependence of the measured galaxy 
power spectrum in redshift space. We highlight that, in the 
presence of either nonlinearity or a stochastic relation, the 
correspondence between the coefficients of /x 2 and fj 4 , and 
/lt is no longer so simple and a more involved treatment is 
required. 
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